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I. INTRODUCTION 

The asymptotic eigenvalue distribution of partial dif- 
ferential equations such as the Schrodinger equation for 
free particles or the one-dimensional Helmholtz equation 
for sound waves has been of interest as early as from 1912 
on when Hermann Weyl and Richard Courant first stud- 
ied the problem [TJ [5] . They found expressions for the 
asymptotic level number N(k) in closed systems to be 
proportional to k d , with d the spatial dimension of the 
system. The so called "Weyl law" , which has been well 
known from then on [3J, states that for closed quantum 
systems, every accessible Planck cell in phase space is oc- 
cupied by one quantum state. A generalization to chaotic 
open systems, where complex resonances k n = k n — T n /2 
with mean energies k n and lifetimes T„ replace real eigen- 
values k, has been proposed in the 1990s [H [5]. The 
number of resonances 

NQt) = {k n : Re(k n ) < k; Im(fc„) > -C} (1) 

inside a rectangle in the complex plane defined by the 
energy k and the strip width C is conjectured to be pro- 
portional to k a with the exponent 



D + l 



(2) 



being related to the non-integer fractal dimension D of 
the chaotic repeller. The number a takes the role of 
the effective number of degrees of freedom. The repeller 
is the set of all classical trajectories that stay trapped 
for t — > +oo or t ^ —oo. Considering only the stable 
manifold W s of trapped trajectories for t — > +oo in a 
suitable Poincare surface of section, the fractal Weyl law 
reads [3] 



N(k) cx k dn - 



(3) 



with g?h the Hausdorff dimension of the chaotic repeller's 
stable manifold. 

The fractal Weyl law ^ has been investigated for var- 
ious two-dimensional systems, e.g. a triple Gaussian po- 
tential [7] , the three-disk billiard [6] , an optical microsta- 
dium resonator [B] and a modified Henon-Heiles potential 




FIG. 1. The four-sphere billiard. Four spheres of equal radius 
R are located on the vertices of an equilateral tetrahedron 
(indicated by bars) with edge length d. Shown is the case of 
d/R = 6. 



[9] as well as for quantum maps, e.g. the kicked rotator 
[T01 ITT] . The systems under consideration so far have all 
been at most two-dimensional. We provide a first inves- 
tigation of a three-dimensional system in this paper. 

The problem under consideration is scattering in the 
four-sphere billiard. This system is characterized by four 
spheres of the same radius R located on the vertices of an 
equilateral tetrahedron with edge length d as visualized in 
Fig. [T] The relevant configuration parameter is the ratio 
d/R. Different from e.g. the two-dimensional three-disk 
billiard, the four-sphere billiard is an open system even 
for the case of touching spheres at d/R = 2. 

The paper is organized as follows. The fractal dimen- 
sion of the repeller is computed in Sec. [n] Periodic-orbit 
theory and cycle expansion methods are applied in Sec. 
|III| to obtain semiclassical resonance spectra. The clas- 
sical escape rate is determined in Sec. |IV| Results which 
allow for a test of the fractal Weyl law are presented in 
Sec. [V] Concluding remarks are given in Sec. VI 



II. GAUGING THE REPELLER 

As the fractal dimension D of the chaotic repeller en- 
ters into the fractal Weyl law ([3]), it is crucial to deter- 
mine D accurately. The fractal structures in the four- 
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sphere billiard have already been studied experimentally 
and theoretically [T2Tll5| . however, all investigations so 
far have been limited to small values of the configuration 
parameter d/R < 2.5. In this paper the fractal dimen- 
sion of the repeller is determined accurately for a wide 
range of the parameter d/R. 



A. Fractal repeller 

The stable manifold W s of the chaotic repeller is a 
fractal in phase space. The time spent in the scatter- 
ing system can be measured by the time-delay function 
T counting the number of reflections experienced by a 
trajectory. Choosing initial conditions in a plane parallel 
to the plane spanned by three of the spheres, it is possi- 
ble to iterate trajectories entering the scattering system 
such that the fourth sphere is visited first. For large ra- 
tios d/R, the boundary to the region that contains those 
initial conditions is the projection of the fourth sphere 
onto the plane from which the trajectories are iterated. 
As there are three distinct possibilities to visit the next 
sphere, there are three regions of higher values of T in- 
side this circle. Repeating this line of argument for any 
region belonging to a given visitation sequence, the struc- 
ture of the fractal repeller can be understood. Figure [2] 
illustrates the fractal structure. 

In principle, the Hausdorff dimension dfj can be calcu- 
lated from box-counting. This procedure, however, is not 
suited for the billiards under consideration as it requires 
iteration of a vast amount of initial conditions on a grid. 
The regions of high values of T that exhibit fractal prop- 
erties may not be resolved with acceptable computational 
effort. A method better suited to billiards is introduced 
below. 



B. Estimating du through Hausdorff sums 

Finite numerical precision and finite computing time 
available prevent the determination of initial conditions 
that lead to trapped orbits. Instead, it is possible to es- 
timate the Hausdorff dimension from regions of finite T. 
We introduce the auxiliary quantity An that denotes the 
area of the n-th region of initial conditions with T > i 
reflections in the surface of section, and define the quan- 
tities 

A-«( S ):=]T(aW) S (4) 

n 

which will be called Hausdorff sums below. These sums 
have the following properties [HI [17] : 

!oo for < s < du 

const. > for s = rf H (5) 
for du < s < oo 




-0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 



x 




-0.02-0.015-0.01 -0.005 0.005 0.01 0.015 0.02 



x 

FIG. 2. (a) Time-delay functions in the range 2 < T < 5 
for d/R = 3 and (b) a magnification thereof in the range 2 < 
T < 7. The time-delay functions are drawn as functions of the 
coordinates x, y in the surface of section. The colors indicate 
the value of T and for clear identification some regions are 
explicitly labeled with T. Only the boundaries of individual 
regions are drawn; initial conditions chosen within the regions 
experience the same number of reflections. The self-similarity 
suggests that the stable manifold is a fractal set. 



The property for s = du stems from the fact that the 
Hausdorff sums by definition are smooth functions of the 
variable s. This allows one to estimate the Hausdorff 
dimension du by intersecting K^'(s) for different i [T4l 

Eng. 

Existing methods of estimating the repeller's fractal di- 
mension have been confined to a narrow range of the pa- 
rameter d/R, in particular to the case of almost-touching 
spheres [HI HI] ■ The algorithm discussed in the following 
allows for calculations in a wider range of d/R. To esti- 
mate the stable manifold W s in the surface of initial con- 
ditions, it will be necessary to accurately find boundaries 
of regions of a given visitation sequence of the spheres 
and reflection count. Once both a point inside the region 
in question and a point outside are known, interval bi- 
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section may be used to compute the boundary point on 
the line connecting both points. The bisection condition 
uses the visitation sequence and the reflection count T, 
i.e. the length of the symbolic code. 



1. Finding regions 



The algorithm used in this paper relies on the structure 



of the time-delay functions discussed in Sec. II A In the 
calculations, the following assumptions are made: 

• Regions of a specific order of visits with the scat- 
terers described by the symbolic code are non- 
overlapping. 

• Within a region of a given order of visits, there are 
exactly three more regions each corresponding to 
additional visits at one of the three other spheres. 

Both conditions may be violated for d/R close to 2, i.e. 
the case of almost touching spheres. We find the assump- 
tions to be fulfilled for regions with T > 2 and configu- 
rations d/R> 2.5. 

All steps of the procedure are based on the Poincare 
surface of section chosen such that iteration starts from 
a plane parallel to the plane spanned by the three closest 
spheres. The velocities are chosen parallel to the 2-axis 
such that the uppermost sphere is visited first. Let us 
assume that regions with T m i n < T < T max are sufficient 
for an estimation of dn- Under this assumption, it is 
possible to find regions approximating the repeller with 
the following procedure. 

In the a first step, the projection of the sphere visited 
first onto the surface of section is determined. For small 
d/R, this region is a circle, for larger d/R, projections 
of the other spheres may be cut out of the circle. This 
is done by randomly choosing a point in the surface of 
section and an interval bisection between this point and 
points equally distributed on a large circle fully contain- 
ing the projection of sphere 4. By assumption, inside 
this region there are three other regions with T = 2. The 
corresponding visitation sequences differ in the second 
character. Once a point inside each of the regions with 
T = 2 is found, in a second step, polygonal chains form- 
ing boundaries to each of the new distinct regions are 
calculated. This procedure is iterated until all desired 
regions corresponding to T m ; n < T < T max have been 
found. 



2. Areas from polygonal chains 

One possibility to store the boundaries is by keeping 
a polygonal chain. As in this procedure the number 
•^regions of regions grows exponentially with N reg i ons — 
3 T_1 , this way of data storage is memory-expensive. 
However, the areas enclosed by the polygons are easily 



TABLE I. Numerical values of the Hausdorff dimensions dH of 
the stable manifold W a for various configuration parameters 
d/R. All decimal digits are significant. 



d/R 
dn 



2.5 3 4 5 6 8 10 
0.4774 0.3818 0.2992 0.2596 0.2354 0.2063 0.1888 



calculated using numerical quadrature of the area given 
by 



2tt 



A 



r 2 (ip)dip, 



(6) 



with r(tp) the distance of the boundary point from the 
"midpoint" of the region. A fairly low number of sup- 
porting points has proved to be sufficient for very high 
precision. All calculations have been performed with 101 
supporting points. 



3. Areas from ellipses 

An alternative to the memory-expensive storage of 
polygonal chains is to approximate the boundary by an 
ellipse described by the polynomial 



a-iy 



a^xy + a^x + a^y = 1 



(7) 



For five or more known points {xi,y{) of the boundary 
the coefficients a\, . . . , 05 can be determined from a lin- 
ear least-squares fit. The semi-major axes a, b, the center 
shift (xq, yo) as well as the rotation angle (p of the ellipses 
can be easily extracted from the coefficients of the poly- 
nomial in Eq. ([7]). 

Fitted ellipses allow to improve accuracy as it is now 
possible to shift the regions' "midpoints" used in the con- 
struction of the polygonal chains to the midpoint of the 
ellipses. All bisections for the polygonal chain boundaries 
are repeated in such a way that all lines connecting the 
ellipse's midpoints and the boundary points intersect at 
identical angles. This will be beneficial for the quadra- 
ture of the areas entering into the Hausdorff sums. 

Once all desired boundaries have been calculated, the 
fractal dimension dn can be estimated. To build the 
Hausdorff sums Q, the areas enclosed in the individ- 
ual regions have to be known. From the ellipses fitted to 
the boundaries, the area A is trivially given by 



A = irab , 



(8) 



where a and b are the semi-major axes. 

Calculations have been performed for d/R — 2.5 to 
d/R = 10. Sample plots for intersected Hausdorff sums 
(s) are shown in Fig. [3] Results for the Hausdorff 
dimension dn are compiled in Fig. [4] and in Table [I] The 
calculations using polygonal chains agree up to four dec- 
imal digits with the calculations using fitted ellipses. 
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FIG. 3. Intersected Hausdorff sums K (i \s) for various re- 
flection numbers i — T calculated from (a) a polygonal chain 
with 101 supporting points and (b) ellipses fitted to polygonal 
chains for d/R = 6. As can clearly be seen, the intersection 
points for all shown curves agree perfectly. For this reason, 
the Hausdorff dimensions can be determined to a precision of 
at least 4 significant digits. 



FIG. 4. Hausdorff dimension du of the stable manifold W s as 
function of the ratio d/R. All data points have been obtained 
by intersecting Hausdorff sums as demonstrated in Fig. [3] 



niques of semiclassical quantization presented below are 
better suited for billiard systems. 

As one of the great achievements of semiclassical 
physics, Gutzwiller's trace formula provides a mean of 
quantizing a system via periodic orbits 21j. Unfortu- 
nately, the trace formula is plagued by serious conver- 
gence problems. In chaotic systems, the number of pe- 
riodic orbits typically grows exponentially with length I, 
and this growth usually cannot be compensated by the 
decrease of the amplitude factors. A method of improv- 
ing the convergence ideally suited for billiard systems is 
based on the Gutzwiller-Voros zeta function. 

The logarithmic derivative of the function 



Z(k) = JJ(fc - k n ) 



(9) 



Figure [4] clearly shows that with decreasing d/R the 
intersection of the stable manifold W s with the Poincare 
surface of section fills the plane denser. The repeller's di- 
mension g?h thus increases as the tetrahedron gets packed 
more densely. 

In summary, the method presented here establishes 
a fast and very precise method of gauging the repeller. 
Though the assumptions are quite strong, they hold over 
a wide range of the ratio d/R. 



III. SEMICLASSICAL RESONANCES 

Studying a billiard system in a purely quantum me- 
chanical fashion turns out to be intricate since - although 
particles move freely in between the scatterers - it is a 
demanding task to find wavefunctions that vanish on all 
scatterer's boundaries simultaneously. For attempts on 
A^-sphere scattering systems in three dimensions, see [15] : 
for the two-dimensional three-disk scattering system see 
[T§] . Results and comparisons of methods for the four- 
sphere scattering system are presented in [50] • The tech- 



with the quantized wavenumbers k n of a billiard system 
yields the density of states 

g{k) = --lm-^-lnZ(fc) = --ImV- } . 

7T ak 7T * — ' k — k n + le 

n 

(10) 

Voros 22J proposed a semiclassical formulation of Z(k) 
which reads for billiard systems 



•^Gv(fc) = exp 



O ^ 



(-1) 



rn p giripfe 



V|det(M/-l)| 



(11) 



where l p is the length of a primitive periodic orbit and n p 
the number of reflections on hard wall boundaries. The 
index r counts the number of repetitions of a primitive 
periodic orbit. The reduced monodromy matrix M p pro- 
vides information on the linear evolution of a small devi- 
ation from an initial condition belonging to a periodic or- 
bit over one period [23] • The eigenvalues of M p quantify 
the stability of the periodic orbit. Due to the symplec- 
tic structure of Hamiltonian mechanics, the eigenvalues 
come in tuples A, 1/A, A*, and 1/A*. 
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TABLE II. Character table for the group T d [20]. The group 
has five irreducible representations: the one-dimensional rep- 
resentations Ai and A%, the two-dimensional representation 
E and the two three-dimensional representations T\ and T 2 . 



T d 


E 8C3 3C 2 6S4 60-d 


A 1 


11111 


A 2 


1 1 1-1-1 


E 


2-1 2 


Ti 


3 0-1 1-1 


T 2 


3 0-1-1 1 



For systems with symbolic dynamics such as bil- 
liards, the method of cycle expansion |24H27j has proved 
to be especially successful. A cycle expansion of the 
Gutzwiller-Voros zeta function Zcv(k) in Eq. (Ill is 
achieved by replacing (— l) rn p in Eq. ( |lTj ) by the term 
(— z) rnp depending on the book-keeping variable z, ex- 
panding Zqy as a power series in z and then truncating 
the series. The highest power of z equals the maximum 
cycle length n max contributing to the cycle expansion. 
After truncation, z has to be set to z = 1. The cycle- 
expanded zeta function has better convergence behavior 
over the trace formulas as individual terms tend to can- 
cel. 



A. The symmetry group Td 

The four-sphere billiard has discrete tetragonal sym- 
metry [20] . The associated symmetry group T& contains 
all symmetry operations that leave a regular tetrahedron 
invariant. In particular, there are the identity operation 
E, 4 rotations C3 by 27r/3 around the axes defined by a 
vertex of the tetrahedron and the center of the facing tri- 
angular boundary surface, 4 more rotations C3 2 by 4-7r/3 
around the same axes, 3 rotations C2 by it around the 
axes intersecting the middle points of opposing edges, 
6 reflections ad at planes perpendicular to the tetrahe- 
dron's edges and also containing another vertex; and, 
furthermore, 3 permutations of the vertices £4 which can 
be written as a combination of a rotation C4 by 7r/2 and 
a reflection ah at the plane perpendicular to the main 
rotation axis, i.e. the axes of C3. Finally, the symmetry 
group Td also contains 3 distinct three-times repeated 
rotary reflections S4 3 . 

The character table of the symmetry group is given 
in Table [TTJ The symmetry group can be decomposed 
into 5 invariant subspaces, i.e. the representation matri- 
ces D of Td can be decomposed into block-diagonal form 
where the diagonal elements contain matrices represent- 
ing the group's elements. The representation is called 
"irreducible" if no further decomposition is possible. 

From the character table, statements on the wavefunc- 
tions ip can be made. For the one-dimensional representa- 
tions Ai and A 2 , the effect of the symmetry transforma- 
tions is described by a multiplication with the character 



X- For the totally symmetric A\ subspace, all characters 
are equal to 1, therefore the wavefunctions have the full 
symmetry of T^, i.e. ip is not altered by any symmetry 
transformation. In the A2 subspace, the wavefunction 
changes sign under the reflection ad and permutation S4. 
For representations of higher dimension, the effect of the 
symmetry transformations cannot be described in such a 
simple way. 

We note that repeated application of symmetry trans- 
forms may be identical with other elements of the sym- 
metry group, e.g., ad 2 — C2 2 = C3 3 = S4 4 = E and 
S4 2 = C%. These identities will be useful for the symme- 
try decomposition of zeta functions discussed below. 



B. Symbolic dynamics and periodic orbits in the 
four-sphere scattering system 

The method of cycle expansion requires all periodic 
orbits up to a given maximum cycle length n r . In special 
billiards it is convenient to assign a symbolic code to each 
periodic orbit. For the four-sphere scattering system the 
symbolic dynamics and periodic orbits have already been 
introduced in Ref. [2DJ . For the convenience of the reader 
we briefly recapitulate the central ideas and properties of 
the orbits. 

In the four-sphere scattering system, periodic orbits 
are determined by the periodic sequence in which the four 
scatterers are visited. For most cases, any sequence cor- 
responds to one periodic orbit. If this one-to-one corre- 
spondence is not given, i.e. if some orbits become unphys- 
ical because they penetrate the scatterers, one speaks of 
pruning. This is the case for small center-to-center sep- 
arations. In the four-sphere scattering system, the sym- 
bolic dynamics has been shown to be pruned for config- 
urations with d/R < 2.0482 [20] , 



1. Periodic orbits in the fundamental domain 

In systems with discrete symmetries such as the four- 
sphere billiard, which is invariant under all symmetry 
operations of the tetrahedron group Td, whole classes of 
orbits are equivalent to each other, e.g., the 6 orbits which 
are scattered back and forth between two spheres can be 
mapped onto each other using symmetry operations of 
Td- Furthermore, cyclic permutation of the sequence of 
spheres leaves the orbits invariant. For these reasons, it 
is appropriate to use the symmetry properties to intro- 
duce the following short notation [2D]. First, define the 
plane of reflection as the plane that contains the centers 
of the last three distinct spheres visited. Then, instead of 
labeling all spheres individually, the label will be used if 
the orbit visits the last sphere once more, 1 will indicate 
a visit at the third other sphere in the same plane of re- 
flection, whereas the label 2 will be used for a visit at the 
fourth sphere outside the plane of reflection. With this 
nomenclature, all symbolic codes containing the charac- 
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TABLE III. Primitive periodic orbits up to cycle length np = 
2 for d/R = 4. The reduced symbolic code p as well as the 
symmetry hp of each cycle is given. Furthermore, the real and 
imaginary parts of the stability eigenvalues A-*' are tabulated. 
All numbers have been rounded to five decimal digits. The 
shortest cycle, labeled by 0, which visits two spheres in turns 
has ambiguous symmetry. Both the rotation about n, C2, as 
well as the reflection about the plane perpendicular to the line 
connecting the sphere's center, ad, map this particular orbit 
onto itself. 



V 



Re A 



(i) 



ImA 



(i) 



Re A. 



(2) 



ImA 



(2) 



<7 d ,C 2 2.00000 5.82843 0.00000 5.82843 0.00000 

1 C 3 2.26795 -7.09669 0.00000 5.75443 0.00000 

2 S 4 2.31059 -3.11111 5.69825 -3.11111 -5.69825 

01 a d 4.34722 -46.21054 0.00000 32.08725 0.00000 

02 C 3 4.35831 -14.95013 35.68205 -14.95013 -35.68205 
12 S 4 4.58593 43.79192 0.00000 -39.51750 0.00000 



ter 2 are three-dimensional, whereas orbits corresponding 
to sequences of and 1 are two-dimensional. The new 
labeling reduces the number of characters in the alphabet 
to three, i.e. the code is ternary. This reduction corre- 
sponds to a reduction of the full physical phase space 
M to the fundamental domain M from which the whole 
phase space can be reconstructed by applying the symme- 
try group's elements. Note that the symmetry reduced 
orbits are in general shorter than the corresponding phys- 
ical ones. Only the symmetry reduced orbits that have 
the identity operation E as maximum symmetry have the 
same length as the corresponding physical orbits. The re- 
duced orbits of symmetry classes ad and C2 yield twice 
as long physical orbits, C3 orbits are three times longer 
than in the physical space, and, finally, S4 orbits have 
quadruple length. 



2. Finding periodic orbits 

Periodic orbits of the four-sphere scattering system are 
calculated by varying, for a given symbolic code, the re- 
flection points on the spheres until the length of the or- 
bit takes its minimum value. Details of the numerical 
periodic orbit search are described in [30] and the com- 
putation of the monodromy matrix is explained in Refs. 
[2"51 l2l?] . Table III lists the first few periodic orbits in the 



fundamental domain as well as their properties for the 
ratio d/R — 4. The table also gives the maximum sym- 
metry operation that leaves the orbit invariant. In the 
fundamental domain, this operation corresponds to the 
operation that maps the endpoint of the orbit in the fun- 
damental domain onto the starting point. Note that, for 
example, the cycle io which visits two spheres in turns 
is periodic in the fundamental domain, but not in full 
physical space. In the full domain, the start and end- 
point of the 0-cycle are not identical; application of the 
rotation C2 respectively the reflection ad yields back the 



full periodic orbit. 



C. Discrete symmetries and cycle expansion 

In systems with discrete symmetries the full physical 
spectrum can be decomposed into spectra belonging to 
different representations of the symmetry group. The 
discrete symmetries lead to symmetry factorized zeta 
functions, which allow for the computation of quantum 
spectra belonging to a specific symmetry subspacc. The 
symmetry decomposition of zeta functions has been elab- 
orated by Cvitanovic and Eckhardt (30] and examples 
have been given for the symmetry groups of various two- 
dimensional A^-disk pinball models. Here, we present ex- 
plicit results for the tetrahedron group Td of the three- 
dimensional four-sphere scattering system. As in [30] we 
first discuss the symmetry decomposition of the dynam- 
ical zeta function 



Z{k)=\{(l-t p (k)) 



(12) 



which is obtained from the Gutzwiller-Voros zeta func- 



tion (11) with the approximation 



|det(M/ - l)f 1/2 » \\^\^\- r/2 = e -K 1) +^ 2) )-/2 

(13) 



and the definition 

t p (k) = e i^-« ) +4 2) )/2 



(14) 



and then generalize the results for the symmetry decom- 
position of the Gutzwiller-Voros zeta function. 

In quantum mechanics, the full Hilbert space H of the 
problem factorizes into subspaces belonging to certain 
irreducible representations of the symmetry group, i.e. 



(15) 



for the four-sphere scattering system. In |30j it is pointed 
out that zeta functions can be factorized in a similar way. 
The fundamental domain of phase space is sufficient for 
all computations, as the whole phase space M can be 
obtained from the fundamental domain M by 



M 



(16) 



heG 



where G is the symmetry group. Evaluating traces of 
transfer operators in the fundamental domain M, this 
symmetry reduction results in [26} 130] 

{l-t p ) m * = det (1 - £> (ftp) tp) , (17) 

with Trip the multiplicity of a primitive cycle p. 

These expressions could be evaluated using a certain 
explicit representation D a (h) of the group's symmetry 
operations h. However, this is a computationally rather 
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demanding endeavor. Instead, the determinants can be 
expressed in terms of traces \ which can be read off from 
the symmetry group's character table (see Table [TTJ) . For 
example, the expansion of det(l — D(h)t) for dimension 
d = 3 reads 

det(l - D(h)t) = 1 - x (h)t + \ ( X (hf - x (h 2 )) t 2 

+ l(x(hf-3x(h)x(h 2 )+2 X (h 3 ))t< 

(18) 

where the trace of D(h) is as usual denoted by x(h)- 
Carrying out this procedure explicitly, one obtains the 
factorizations given in Table |TV} Thus, the zeta function 
in Eq. ( 12 ) can be rewritten in a symmetry reduced ver- 
sion 



Z a =]J(1- D a (hp)t p ) 



(19) 



for the subspace a. The zeta function now depends only 
on the fundamental cycles p. By this procedure, a fac- 
torization 



z(k) = l[z a (k) d ° 



(20) 



If this is possible, a way to use the p for repetitions as 
well has been found. As an example, for the contribution 
of the r-times repeated cycle p to the A\ spectrum, we 
need to solve 



(1 - z r t p . r ) = . 



(23) 



which is true for z r = 1. Thus, in the A\ subspace, all 
weight factors are wp^ r = 1. By this choice, the symmetry 
'factorization is retained. As another example, consider 
the E subspace for cycles with C3 symmetry. Here, so- 
lutions to the equation 







are needed. A factorization is given by 



1 



270773 X 

If. 



1 



2707-/3 x 



0, 



(24) 



(25) 



where the exponentials are the roots z;. Evaluating 
the sum z\ + z 2 r , we find the weight factors wp yr — 

— 1, —1, 2, —1, — 1 . . . for r = 1, 2, A short notation 

for this sequence is given by Wp^ r = 2 cos(27rr/3). By sim- 
ilar calculations, the weight factors wp^ r given in Table |V| 
are determined. 



is achieved. The zeta function Z factorizes into zeta func- 
tions belonging to certain irreducible representations a of 
the symmetry group. The dimensions d a of the represen- 
tations enter into the full zeta function - and with them, 
the quantum multiplicities of resonances belonging to a 
certain subspace. 



1. Assigning weight factors 

The method of cycle expansion expands the zeta func- 
tion Z into a truncated series in which all cycles up to 
a certain cutoff length enter [2"4Tl27| . However, besides 
the primitive cycles, also multiple traversals contribute. 
Therefore, it needs to be clarified how repeated revolu- 
tions can be taken into account. Let us assume the prim- 
itive fundamental cycles p are known. Then, the contri- 
bution of an r-times repeated revolution to the symmetry 
reduced zeta function ( 19 1 is given by polynomials such 
as 



(1 Z £p,r) , 



(21) 



where a dummy variable z has been introduced. The cy- 
cle weights tp >r have the form of the terms in (11) and 



are thus easily calculable from the cycle weight tp of the 
primitive fundamental cycle. By using the factorizations 
given in Table |TV} it is possible to determine the weight 



factor wp tr (k) as the sum of all roots Zi" of the polyno- 
mials given in the table, 



u p,r 



E- 



(22) 



2. Ambiguous symmetry 

The shortest cycle labeled by in the four-sphere sys- 
tem has ambiguous symmetry. It is possible to map this 
cycle onto itself by both the rotation C2 and the reflec- 
tion CTrf. This ambiguity requires special care in the sym- 
metry decomposition. This is particularly important as 
the O-cycle is one of the fundamental cycles that act as 
building block for longer cycles in the sense of cycle ex- 
pansion. The group theoretical weight of the O-cycle can 
be written as [3D] 



ho 



C 2 



(26) 



The symmetry factorization can thus be not one of those 
given in Table |IV| However, it is possible to use a fac- 
torization that contains factors in such a way that the 
factorization is at most the greatest common divisor of 
the factors given for C 2 and aj, in Table [TV] The factor- 
izations and weight factors wo, r are given in Table VI 



With that factorization the product ( 20 1 of all zeta func- 



tions belonging to the irreducible representations of the 
symmetry groups coincides with the cycle expansion ( 12 ) 
using all orbits in the full domain. 

The final form of the Gutzwiller-Voros zeta function 
we use for our calculations is 



Z G v- a (k) = exp 



EE 

p r—l 



1 fflp,i-;a(-z)™' Q lrl i> k 



(27) 



TABLE IV. Symmetry factorization of the zeta function Z for all five irreducible representations of the group Td- The table 
entries give the contribution of each fundamental cycle p to the Euler product Z — Ylp(l — tp). This factorization allows the 
computation of quantum spectra for each symmetry subspace. 
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TABLE V. Weight factors wp )T for r traversals of the primitive 
cycle p. These factors allow for symmetry factorizations with 
repetitions of primitive fundamental cycles. 
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TABLE VI. Factorizations of Z(k) and weight factors wo,r 
for the fundamental cycle with ambiguous symmetry classes 
C 2 , Od in all subspaces a of the symmetry group Td . 
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with p the primitive symmetry reduced cycles, r the num- 
ber of repetitions of it and a the symmetry subspace. 
A symmetry reduced version of the cycle expansion is 
obtained by expanding Eq. (27) into a power series in 



z which is truncated at a maximum cycle length n n 
Then, z has to be set to z — 1. 



D. Harmonic inversion method 



The Gutzwiller-Voros zeta function Z(k) in Eq. (27) 
contains all energy eigenvalues k as complex zeros, and, 
in principle, it is possible to obtain spectra by a numeri- 
cal root search. This method has been successfully used 
for billiards, see e.g. [27]. However, the root search in 
cycle expansions of high order is numerically expensive. 
For statistical purposes it is important not to miss any 
resonances in the strip of the complex plane under consid- 
eration. Therefore, a dense grid of initial root guesses has 
to be used for the root search. Consequently, many reso- 
nances will be found several times. Thus, the problem is 



to distinguish for every new root whether a new distinct 
resonance has been found or if the new zero has already 
been computed. As the number of resonances enters into 
the fractal Weyl law Q through the counting functions 
N(k), it is crucial to count individual resonances only 
once. 

An alternative to the computation of zeros is based on 
the harmonic inversion method for high-resolution spec- 
tral analysis [3lT433j . When Eq. (10 1 is evaluated along 
a line of real- valued k or a shifted line k + iS with real k 
and J, we obtain a spectrum 



9(h) 



E 1 



r„/2 + s 



* (k-k n y + (r n /2 + sy 



(28) 



which is a superposition of resonances with a Lorentzian 
shape. For negative shifts S the Lorentzians are located 
at the positions k n , but with reduced widths T n + 26. 
The basic idea is now to reformulate, via a Fourier trans- 
form, the problem of extracting eigenvalues as a signal 
processing task. Details of the method are given in [S]. 

The procedure of calculating quantum spectra is sum- 
marized as follows: First, the spectrum g(k) is calculated 
as a superposition of Lorentzians. We use the cycle- 
expanded zeta function Z(k) for this purpose. The quan- 
tity 

9(k) = - 1 -^l k lnZ(k) = - 1 -lra^ (29) 

is evaluated along lines parallel to the real axis with dif- 
ferent shifts 5. Thus, the shifts that allow for better 
results in harmonic inversion enter into the cycle expan- 
sion. Then, harmonic inversion is used to obtain the 
Lorentzians' parameters k n and T n for spectra calculated 
with different shifts. In the next step, the spectra are 
filtered via the amplitudes. The quantity g(k) given in 



( 29 ) should give resonances with an amplitude of A n = 1. 
True resonances with amplitudes A n rs 1 can be clearly 
separated from spurious resonances with nearly zero am- 
plitudes. Finally, the spectra for different shifts 6 are 
joined such that the individual strips do not overlap. 



IV. CLASSICAL ESCAPE RATE 

The classical escape rate 70 can be interpreted descrip- 
tively as follows (34] : presume the scattering system un- 
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der consideration is located in a box much larger than 
the system itself. Conducting Nq scattering experiments 
with the same incident energy k, but different incident 
directions, one finds that the number N t of trajectories 
that are inside the box after the time t has passed decays 
exponentially as 



N t oc iV e" 7ot . 



(30) 



The relation of the escape rate and the imaginary part 
of the quantum resonances can be understood from the 
correspondence principle. The number of classical trajec- 
tories inside the box corresponds to the quantum proba- 
bility density (ifi\ip). The decay of this probability, 



-rt 



(31) 



relates to the decay ( 30 ) of the number of classical tra- 



jectories inside the box. Thus, in the classical limit 

Tni/r„ = -— -.■ — (32) 



To 
" 2 



holds. 

The classical escape rate can be calculated by the 
method of cycle expansion as well [S3]. The escape rate 
7o is found to be the largest real zero of a dynamical zeta 
function 



z( S )=i[(i-t P (s)), 



(33) 



with p the primitive periodic cycles and tp the cycle 
weights. For a three-dimensional system, 



tp(s) = 



(34) 



The quantities are the leading stability eigenvalues. 
A generalization to a zeta function for three dimensions 
and multiple traversals r of the primitive cycle p is given 
by 



1 



Z( S )=exp -EE; 



(35) 



with A ~ the leading two stability eigenvalues of p. This 
zeta function can be cycle-expanded as described in Sec. 
|III| Results for the escape rate 70 at various configura- 
tions d/R are given in Table VII 



V. RESULTS 

The fractal Weyl law has been put to test for billiard 
systems before. In 6j, the 3-disk billiard has been stud- 
ied. To make our own results comparable to those given 
in [BJ, we carry out a similiar discussion. 



TABLE VII. Classical escape rates 70"' in order n of the cycle 
expansion for various values of the configuration parameter 
d/R. 

d/R \ 7^ 7^ 7^ 7^ 
4 1.16655 1.16459 1.16440 1.16440 
6 0.85042 0.84977 0.84974 0.84974 
8 0.68259 0.68230 0.68230 0.68230 
10 0.57634 0.57619 0.57619 0.57619 



A. Defining a scale for the strip widths 

For the 3-disk system discussed in [BJ , the strip widths 
C have been chosen in relation to the classical escape 
rate 70. For large values k — ¥ 00, the imaginary part 
of quantum resonances converges to lmk = —70/2 [6, 
I35J. Thus, the discussion of the results is simplified by 
rescaling the strip widths C to 



C :-- 



C 
7o72' 



(36) 



which defines a universal scale independent of the sym- 
metry subspace and the ratio d/R. Similar to [BJ, we 
evaluate the fractal Weyl law for scaled strip widths 
Cg[1;1.6]. 



B. Counting resonances 

We have computed spectra for various values of d/R 
in all symmetry subspaces. Generally, we find the best 
convergence behavior of cycle expansions for large values 
of d/R. Furthermore, the one-dimensional group repre- 
sentations A\ and Ai yield the largest number of con- 
verged resonances. The two-dimensional representation 
E and the three-dimensional T\, T2 representations con- 
verge not as well in cycle expansion since the shadowing 
of individual cycles is less efficient for the weight factors 
of these subspaces. For A\, where all weight factors are 
equal to 1, the best convergence is observed. 

It is important to note that for the tests of the fractal 
Weyl law we have used only converged resonances. For 
example, the A\ resonances at d/R = 10 are sufficiently 
converged in the region Refc < 6000, lmk > —0.45 so 
that counting functions N(k) obtained in orders 10 and 
11 of the cycle expansion fully agree. Thus, comparisons 
with the fractal Weyl law are not affected by the order 
of the cycle expansion. 

The fractal Weyl law ([3| suggests that the counting 
functions N(k) obey a power law, 



N(k) oc k a , 



(37) 



thus, in a logarithmic plot of N(k), straight lines of slope 
a = 1 + dn are expected. A sample spectrum and cor- 
responding counting functions for d/R = 10 are shown 
in Fig. [5j This figure is generic in structure, i.e. we have 
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FIG. 5. (a) Spectrum and (b) counting functions N(k) for 
Ai resonances calculated from cycle expansion in order 11 for 
the ratio d/R — 10. Several counting functions for different 
strips C are shown. The curves can be used as "raw data" 
to fit power laws. In this way, the exponent a in the fractal 
Weyl law can be obtained and compared to the classical cal- 
culations. More than 50 000 resonances have been used in the 
analysis. 



found similar behavior of N(k) in other subspaces and 
for other ratios d/R as well. Thus, a brief discussion of 
these features will be given in the following. 

We first note that the strip width C has to be suffi- 
ciently large, since otherwise, counting would not involve 
reasonably large numbers of resonances. In the spectrum 
shown in Fig. [5ja), we have found converged resonances 
in the relevant strip C e [1; 1.6] for Refc < 6000. For 
small strip widths such as C — 0.28 C = 0.97, it is 
evident from Fig. [BJb) that counting involves only a lim- 
ited number of resonances. Larger strip widths involve 
more resonances in the counting. However, choosing the 
strip width too large, the counting may also involve res- 
onances that may not have been converged. Taking the 
asymptotic behavior of the resonances' imaginary parts 
into account, choosing rescaled strip widths in the inter- 
val C € [1.0; 1.6] turns out to be a reasonable choice. 

Figure [5jh) reveals that the counting functions N(k) 
deviate from power laws that led to straight lines in the 
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FIG. 6. Exponents a obtained from least-squares fits of a 
power law to measured counting functions for (a) Ai reso- 
nances and (b) A2 resonances calculated for d/R = 6 in order 
13 of the cycle expansion. The power law has been fitted to 
the interval k £ [0;2000]. The vertical dotted line gives the 
classical exponent a = 1 + du — 1.2354. 



plot. From this observation one infers that the exponent 
a will clearly depend on the fc-range one fits to. We follow 
[6] and choose the largest interval converged resonances 
have been computed for. 



C. Putting the fractal Weyl law to test 

To provide several tests for the fractal Weyl law, we 
will provide plots showing the exponents a obtained from 
least-squares fitting a power function N(k) cx k a to the 
counting functions calculated from the spectra for vari- 
ous subspaces, fitting ranges [0; fc max ], configuration pa- 
rameters d/R and strip widths C. We have performed 
least-squares fits to match the power function p7| ) to the 
measured N(k). 



1. Configuration d/ R = 6 

For d/R = 6, we obtain the exponents shown in Fig. [6] 
For the Ai subspace, we find a very good agreement for 
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FIG. 7. Exponents a obtained for d/R = 8 from least-squares 
fits of a power law to measured counting functions for (a) A\ 
resonances and (b) A2 resonances. The power law has been 
fitted to several fc-intervals. The vertical dotted line gives the 
classical exponent a = 1 + du — 1.2063. 



FIG. 8. Exponents a obtained for d/R — 10 from least- 
squares fits of a power law to measured counting functions 
for (a) A\ resonances and (b) A2 resonances. Results for sev- 
eral fc-intervals are shown. The vertical dotted line gives the 
classical exponent a = 1 + du — 1.1888. 



moderate values of C < 1.4. The relative error in this 
C-interval is less than 2 percent. However, in the A 2 
subspace, all computed exponents are too large by about 
8 percent for the same C-interval. One possible reason is 
that the fc-range used for fitting is too small. 



2. Configuration d/ R = 8 

Performing the same procedure for a configuration pa- 
rameter of d/R = 8, we obtain the plots shown in Fig. [7] 

Both plots reveal a clear tendency to obey the fractal 
Weyl law within a smaller error range when the range of 
k values used for the fit increases. However, for reasons 
of convergence, longer spectra have not been used. We 
note that for C < 1.4 and k € [0;3000], the error is less 
than 7 percent for the A\ resonances. The exponents 
obtained from A2 resonances are larger than the expected 
exponent. For C — 1.3, the relative error is about 15 
percent. 



3. Configuration d/R — 10 

Finally, the system configuration given by d/R = 10 
has been studied. Results are shown in Fig. [HJ 

The counting functions for the A\ subspace once more 
tend to give the expected exponent when the fc-interval 
used for the fit increases. For k 6 [0;6000], the error is 
less than 3 percent. Again, the A 2 spectra yield expo- 
nents that are too large. Possibly the fc-range investi- 
gated here is not large enough to exhibit the asymptotic 
behavior clearly. 



D. Discussion 

While the classical calculations for the fractal dimen- 
sion du are accurate to at least four significant digits, 
the agreement of the exponents a is at best 2 to 7 per- 
cent for the A\ spectra. However, we note that for the 
two-dimensional three-disk billiard, the errors in the ex- 
ponents are about 5 to 10 percent [6]. Therefore, we 
conclude that the fractal Weyl law for the four-sphere 
scattering system is confirmed with roughly the same ac- 
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curacy as for the three-disk billiard. 

A very large fc-range seems to be necessary for a proper 
investigation. This tendency is also visible for the A 2 
spectra. The exponents obtained from the A 2 spectra 
are too large. However, using larger fc-intervals, the ex- 
ponents seem to approach the correct value for large strip 
widths C. Possibly, if larger spectra were available, the 
expected exponents could be obtained. Unfortunately, 
we are limited by the convergence of the cycle expansions 
we use. The higher-dimensional symmetry subspaces E, 
Ti and T 2 could not be used to put the fractal Weyl to 
test law since the spectra did not contain enough con- 
verged resonances. 



VI. SUMMARY AND OUTLOOK 

This paper provides a test of the fractal Weyl law for 
a three-dimensional scattering system. The four-sphere 
billiard was investigated both classically and quantum 
mechanically. 

In Sec. |nj we have developed a fast and very precise 
method to gauge the repeller. We found estimates for the 
Hausdorff dimension with a relative accuracy of 10~ 4 . 
Although the algorithm is based on strong assumptions, 



it works over a wider range of the configuration parame- 
ter d/R than existing methods. 

In Sec. |III| we have discussed the methods of semiclas- 
sical quantization. We have applied the method of cycle 
expansion to the four-sphere billiard. Furthermore, for 
the first time, the method of symmetry decomposition 
was demonstrated for the Gutzwiller-Voros zeta function 
of the system. 

We have given results in Sec. [V] We have provided 
tests of the fractal Weyl law for various configurations of 
the system. Although we have found the counting func- 
tions N(k) to deviate from power functions, we could 
confirm the fractal Weyl law for the Ai resonances of the 
four-sphere scattering within a small error range. For 
those spectra we did not find a convincing agreement 
of calculated level numbers N(k) with the prediction 
N(k) oc k 1+dm for, there is hope that larger spectra would 
lead to the expected exponent. We also assume that the 
deviations from pure power laws are due to the fact that 
the energy range under consideration is too small. 

As an outlook, the physical origin of the modulations 
in the counting functions N(k) will have to be investi- 
gated. Moreover, it is desirable to study further three- 
dimensional scattering systems to find out to what extent 
the results found for the four-sphere billiard carry over 
and are generic. 
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